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>■ . Abstract 

, We propose a new experimental technique for cyclic voltammetry, 

\^ I based on the first-order reversal curve (FORC) method for analysis of 

■ systems undergoing hysteresis. The advantages of this electrochemi- 
cal FORC (EC-FORC) technique are demonstrated by applying it to 

I dynamical models of electrochemical adsorption. The method can not 

only differentiate between discontinuous and continuous phase tran- 
sitions, but can also quite accurately recover equilibrium behavior 
from dynamic analysis of systems with a continuous phase transition. 
' Experimental data for EC-FORC analysis could easily be obtained 

^ ■ by simple reprogramming of a potentiostat designed for conventional 

Q , cyclic-voltammetry experiments. 
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1 Introduction 



Recent technological developments in electrochemical deposition have made 
possible experimental studies of atomic-scale dynamics [1]. It is therefore 
now both timely and important to develop new methods for computational 
analysis of experimental adsorption dynamics. In this paper we apply one 
such analysis technique, the first-order reversal curve (FORC) method, to 
analyze model electrosorption systems with continuous and discontinuous 
phase transitions. We propose that this electrochemical FORC (EC-FORC) 
method can be a useful new experimental tool in surface electrochemistry. 

The EC-FORC method was originally conceived [5] in connection with 
magnetic hysteresis. It has since been applied to a variety of magnetic sys- 
tems, ranging from magnetic recording media and nanostructures to geo- 
magnetic compounds, undergoing rate-independent (i.e., very slow) magne- 
tization reversal [3]. Recently, there have also been several FORC studies 
of rate- dependent reversal [HEIE]- Here we introduce and apply the FORC 
method in an electrochemical context. 

The FORC analysis is applied to rate-dependent adsorption in two-dimensional 
lattice-gas models of electrochemical deposition. We study the dynamics of 
two specific models, using kinetic Monte Carlo (KMC) simulations. First, 
we consider a lattice-gas model with attractive nearest-neighbor interactions 
(a simple model of underpotential deposition, UPD), being driven across 
its discontinuous phase transition by a time-varying electrochemical poten- 
tial. Second, we study a lattice-gas model with repulsive lateral interactions 
and nearest-neighbor exclusion (similar to the model of halide adsorption on 
Ag(lOO), described in Refs. [11 [HI El 110]), being similarly driven through its 
continuous phase transition. Some preliminary results of this work have been 
submitted for publication elsewhere [llj. 

The rest of this paper is organized as follows. In Sec. [2]the FORC method 
is explained. The lattice-gas model used both for systems with continuous 
and discontinuous transitions is introduced in Sec. [31 In Sec. [H the dynamics 
of systems with a discontinuous phase transition are studied. The dynamics 
of systems with a continuous phase transition are studied in Sec. [3 Finally, a 
comparison between the two kinds of phase transitions is presented together 
with our conclusions in Sec. [61 
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Figure 1: (a) Schematic diagram of two first-order reversal curves (FORCs) 
9{fir,fti), separated by the reversal-field step size Aflr (solid lines). The 
dotted lines represent the major hysteresis loop, (b) The corresponding 
voltammetric currents, i{fir,fJ'i) oc dO/djii. 

2 The EC-FORC Method 

For an electrochemical adsorption system, the FORC method consists of 
saturating the adsorbate coverage ^ in a strong positive electrochemical po- 
tential /i (proportional to the electrode potential, with the same sign for 
anions and opposite sign for cations) and, in each case starting from satu- 
ration, decreasing /i at a constant scan rate to a series of progressively 
more negative "reversal potentials" Jir- (See Fig. [1]). Subsequently, the po- 
tential is increased back to the saturating value at the same rate VL 0. The 
method is thus a simple generalization of the standard cyclic voltammetry 
(CV) method, in which the negative return potential is decreased for each 
cycle. This produces a family of FORCs, O^firyfii), where 9 is the adsorbate 
coverage and fii is the instantaneous potential during the increase back to- 
ward saturation. Alternatively, as it is usually done in CV experiments, one 
can record the family of voltammetric currents, 

</x.,/x.) = -7e^^, (1) 

where 7 is the electrosorption valency [121 [131 E] and e is the elementary 
charge. Although we shall not discuss this further here, it is of course also 
possible to fix the negative limiting potential and change the positive return 
potential from cycle to cycle. 
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Figure 2: A subset of the family of FORCs. The circled points are an example 
of the 'slanted' 5 x 11 array of data points used to calculate p from 6{jir, Pi)- 
The filled squares represent a 'square' 5x5 array and the hollow squares 
represent a 'square' 5 x 11 array. See discussion in the text. 



The family of FORCs consists of N FORCs measured at values of jlr that 
are evenly spaced by A/i^. Each of these FORCs has data points measured 
at the (in general different) constant spacing A/2j [3]. It is further useful to 
calculate the FORC distribution, 

2djirdjii 27e(d/ij/dt) djir 

which measures the sensitivity of the dynamics to the progress of reversal 
along the outermost hysteresis loop or major loopQ 

The details of the calculation of p depend on whether the measured 
quantity is the coverage 6{pr,Pi), which is the most convenient in simu- 
lation studies, or the voltammetric current i{fir,Pi), which would be most 
usual in CV experiments. In the latter case, the calculation involves sim- 
ple one-dimensional numerical differentiation with respect to pr, e.g. p = 

e{djli / dt)]~^[i{flr + Aflr,Pi) — i{pr, Pi)]/ Apr- Geometrically it is propor- 
tional to the vertical distance between the current traces. (See Fig. [Hb).) 

^Note that to normalize the FORC distribution, the term ^S{p,i — Mr ) ^^^^^ I p . 
must be added to Eq. (U) [Kj. Here we consider the distribution only away from the line 
fii — fir- The additional term could be found from the major loop. 
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In the case that O^firyfii) is the measured quantity, as in the simulations 
presented here, the calculation of p is somewhat more complicated. An array 
of data points 9{fir,Pi) from consecutive FORCs is used to calculate p as 
shown in Fig. [2l A polynomial surface is fit to this array of data points to 
provide a functional form of 9{pr, Pi), the mixed partial derivative of which is 
proportional to p as given by Eq. ([2]). Away from the major hysteresis loop, 
the array of data points used to calculate p can be taken to be 'square', i.e., to 
start and end at the same value of pi for each of the FORCs involved (squares 
in Fig. [2]). Close to the major loop, this is not geometrically possible. As 
shown in Fig. [2], we therefore used a 'slanted' array of data points to calculate 
p. To obtain the values of p at grid points other than the center point of 
the grid {pi, closer to pr) we used a third-order polynomial to fit 9{pr,Pi). 
This allows the evaluation of the mixed partial derivative p at all points 
of the grid, while a second-order polynomial fit only allows the evaluation 
of p at the center point. For consistency and convenience, a third-order 
polynomial was used with the three kinds of grids described above and at 
all points close to and away from the major loop. The three kinds of grids 
were compared for points away from the major hysteresis loop. While the 
larger 5x11 grids ('square', and 'slanted') produced very similar and smooth 
FORC diagrams, the FORC diagram produced by the 5x5 grid was more 
noisy. Consequently, in the results shown in this paper, the 5 x 11 'slanted' 
grid was used to calculate p from 9{pr,Pi) for all values of the arguments. 

Another approach to calculating p near the major loop would be to use an 
'extended FORC diagram' [15], for which a square grid of points can be used 
even near the major loop, but where the delta function normalization term 
(See Footnote 1) will automatically be included by the polynomial fitting 
procedure described above. The methods for calculating p that are described 
here, are quite straightforward, but also rather slow. A faster method is 
based on the Savitzky-Golay smoothing algorithm [161 [T7]. It is described in 
Ref. [IB]. 

The FORC distribution is usually displayed as a contour plot called a 
'FORC diagram.' A positive value of p indicates that the corresponding 
FORCs are converging with increasing pi, while a negative value indicates 
divergence. The physical significance of the sign of p will become clearer 
when we apply the EC-FORC method to discontinuous and continuous phase 
transitions in Sees. H] and [5], respectively. 
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3 Model 



KMC simulations of lattice-gas models, where a Monte Carlo step (MCS) 
corresponds to one attempt to cross a free-energy barrier, have been used 
to simulate the kinetics of electrochemical adsorption in systems with both 
discontinuous [9l [IHl |19l EO] and continuous [3, [21] phase transitions. The 
energy associated with a lattice-gas configuration is described by the grand- 
canonical effective Hamiltonian for an L x L square system of adsorption 
sites, 

l2 

= - ^ (pijCiCj - fi'^Ci , (3) 

i<j 1=1 

where J2i<j is a sum over all pairs of sites, (pij are the lateral interaction 
energies between particles on the zth and jth sites measured in meV/pair, 
and yU is the electrochemical potential measured in meV/atom. The local 
occupation variables Cj can take the values 1 or 0, depending on whether site 
i is occupied by an ion (1) or solvated (0). The sign convention is chosen 
such that fi > favors adsorption. Negative values of (f)ij denote repulsion, 
while positive values of (pij denote attraction between adsorbate particles on 
the surface. In addition to adsorption/desorption steps, we include diffusion 
steps that have a free-energy barrier comparable to the adsorption/desorption 
free-energy barrier [7|. 

The dynamics of the model are studied by a KMC simulation with the 
computational time unit of one Monte Carlo Step per Site (1 MCSS = 
MCS). For a discussion of the relation between this simulated time unit and 
real, physical time, see Ref. [7]. In each MCS of the simulation, an adsorption 
site is chosen at random and a move (adsorption, desorption, or diffusion) is 
attempted. The transition rates from the present configuration to the set of 
new possible configurations are calculated. A weighted list of the probabilities 
for accepting each of these moves during one MCS is constructed using Eq. ([5]) 
below, and used to calculate the probabilities -R(F|I) of the individual moves 
between the initial state I and final state F. The probability for the system to 
remain in the initial configuration at the end of the time step is consequently 
R{l\l) = l-SF^ii?(F|I) [3 [8]. 

Using a thermally activated, stochastic barrier-hopping picture, the en- 
ergy of the transition state for a microscopic change from an initial state I 
to a final state F is approximated by the symmetric Butler- Volmer formula 
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fT.=^ + A.. (4) 

Here Uj and f/p are the energies of the initial and final states as given by 
Eq. ([3]), Tx is the transition state for process A, and A a is a "bare" bar- 
rier associated with process A. This process can here be either nearest- 
neighbor diffusion (Ann), next-nearest-neighbor diffusion (Annn), or adsorp- 
tion/desorption (Aa/d)- The probability per unit time for a particle to make 
a transition from state I to state F is approximated by the one-step Arrhenius 
rate [221 [23l El] 

- .e.p {-"-I^) ^ .e.p (-^) exp {-"-01) , (5, 

where is Boltzmann's constant and T is the absolute temperature. Here, u 
is the attempt frequency, which sets the overall timescale for the simulation. 
We set u equal to 1 MCSS~^, so that the transition probabilities in a single 
time step (MCS) are given by R{F\l) = 7^(F|I)L-2 MCSS. The electrochem- 
ical potential ft is changed each MCSS, preventing the system from reaching 
equilibrium at the instantaneous value of /2. 

Independent of the diffusional degree of freedom, attractive interactions 
{(pij > 0) produce a discontinuous phase transition between a low-coverage 
phase at low p,, and a high-coverage phase at high fl. In contrast, repulsive 
interactions < 0) produce a continuous phase transition between a low- 
coverage disordered phase for low /i, and a high-coverage, ordered phase for 
high fl. Examples of systems with a discontinuous phase transition include 
underpotential deposition [T^ [201 [SS]; while the adsorption of halides on 
Ag(lOO) [3 [HI [2Sl [2Z] are examples of systems with a continuous phase 
transition. 

4 Discontinuous Phase Transition 

A two-dimensional lattice gas with attractive adsorbate-adsorbate lateral in- 
teractions that cause a discontinuous phase transition is a simple model of 
electrochemical underpotential deposition [191 l20l |25l [28]. Using a lattice- 
gas model with attractive interactions on an L x L lattice with L = 128, 
a family of FORCs were simulated, averaging over ten realizations for each 
FORC at room temperature. The lateral interaction energy (restricted to 
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nearest-neighbor) was taken to be (pij = 0nn = 55meV, where the positive 
value indicates nearest-neighbor attraction. For this value of 0nn, room tem- 
perature corresponds to T = O.STc, where Tc is the critical temperature. 
The barriers for adsorption/desorption and diffusion (nearest-neighbor only) 
were Aa/d = Ann = 150 meV, corresponding to relatively slow diffusion [20] . 
Simulation runs with faster diffusion (Ann = 125 meV) and the same adsorp- 
tion/desorption barrier were little different from the results shown in Fig. [31 
The reversal electrochemical potentials fir associated with the FORCs were 
separated by A/i^ = ImeV increments in the interval [— 200 meV, meV], 
and the field-sweep rate was constant at Q = |A/ij/At| = 0.03 meV/MCSS. 
The FORCs are shown in Fig. [31(a), with a vertical line indicating the position 
of the coexistence value of the electrochemical potential, fio = — llOmeV, 
and a filled circle showing the position of the minimum of each FORC. The 
corresponding voltammetric currents are shown in Fig. [31(b). 

In a simple Avrami's-law analysis, the FORC minima all lie at fii = flo 
[6]. However, in the simulations the minima are displaced. For 6 > 0.5, the 
minima occur at fli < fiQ, precisely at the points where the tendency to phase- 
order, which drives local regions of the system toward the nearby metastable 
state {6 ~ 1), is momentarily balanced by the electrochemical potential, 
which drives the system toward the distant stable state (6* ~ 0). For 9 < 0.5, 
the stable and metastable states are ^ ^ 1 and ~ 0, respectively, and the 
same balancing effect explains the FORC minima occurring at fii > Jiq. The 
net effect is a 'back-bending' of the curve of minima, as seen in Fig. [3]^a). 

In Fig [31(c), the FORC diagram is plotted against the variables Jib = 
{fir + fJ'i)/2 and p,c = (yUr — /ii)/2. These variables are commonly used in 
the literature for plotting FORC diagrams [3], as /2;, denotes the midpoint 
between fir and /ij, and fic is proportional to the distance between these two 
values of /iJl 

The definition in Eq. ([2|) implies that the FORC distribution p should 
be negative in the vicinity of the back-bending. This can be readily seen in 
Fig [31(c). The negative values of p reflect a local divergence of the FORCs 

^In magnetic applications, the variables Hb and He (corresponding to fit and /2c, re- 
spectively) have a clear physical meaning as the bias and coercive fields, respectively, in a 
bistable magnetic system. While this physical meaning does not extend clearly to electro- 
chemical systems, this choice of variables does produce FORC diagrams with less unused 
space than the variables {fir, fii), for which the region fii < fir is forbidden, and so we 
have used it here. However, the analysis described in the following sections could also be 
made with FORC diagrams plotted using the variables {fir, fii). 
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Figure 3: (Color online.) (a) FORCs for a discontinuous phase transition, 
obtained with scan rate = 0.03 meV/MCSS and interactions and barrier 
heights as given in the text. The vertical line shows the coexistence value 
of the electrochemical potential, fx = fio- The minimum of each FORC is 
also shown (filled circles). The thick curve corresponds to the FORC whose 
minimum lies nearest the coexistence value, /ij = /iq. (b) The corresponding 
voltammetric currents, calculated by numerical differentiation of the FORCs. 
See Eq. [H (c) FORC diagram generated from the family of FORCs in (a). 
The positions of the FORC minima are shown as filled circles. The thick, 
straight line corresponds to the FORC marked as a thick curve in (a). 
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Figure 4: The dependence of the FORC minima on the sweep rate fl. The 
figure shows FORC minima for two famihes of FORCs with different sweep 
rates {Q = 0.03 and 0.09 meV/MCSS, respectively). The curves are guides to 
the eye, obtained by smoothing the data using a first-order Savitzky-Golay 
filter with a window of 5 points fl6[ [T7] . 
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away from each other as /ij (and time) increases. This can be considered 
a dynamical instabihty, caused by the competition between the tendency 
to phase-order and the effect of the electrochemical potential. When the 
potential sweep is stopped suddenly at a potential in this unstable region, 
the subsequent time evolution of 9 is non-monotonic: it first approaches its 
metastable value, but then reverses and relaxes reliably to its equilibrium 
value at that potential. The only exception is the point (/I = /2o, d = 0.5) 
along the FORC indicated by bold lines in Fig. [31 It is also interesting 
to note that the curve connecting the minima of the FORCs resembles the 
van der Waals loop in the mean-field isotherm of a fluid system below its 
critical temperature but with an asymmetrical shape about the point 
(/i = /2o, 6 = 0.5) and with a sweep-rate dependent shape as shown in Fig. |H 

5 Continuous Phase Transition 

Using the same Hamiltonian, but with long-range repulsive interactions and 
nearest-neighbor exclusion as appropriate for modeling halide electrosorption 
on Ag(lOO) O El m Ho], KMC simulations were used to produce the family 
of FORCs for a continuous phase transition. The reversal potentials jlr were 
separated by Ajlr = 10 meV increments in the interval [—200 meV, 400 meV]. 
As in Refs. [3 El [ID], the repulsive interactions, with nearest-neighbor 
exclusion and 4>nnn = ~21 meV, are calculated with exact contributions for 
rij < 3, and using a mean- field approximation for r^j > 3. The barriers 
for adsorption/desorption and nearest- and next-nearest-neighbor diffusion 
were approximated on the basis of DFT calculations [2T| as Aa/d = 300 meV, 
Ann = 100 meV, and Annn = 200 meV, respectively [7]. Larger values of 
the diffusion barrier were also used to study the effect of diffusion on the 
dynamics. A continuous phase transition occurs between a disordered state 
at low coverage and an ordered state at high coverage [26l|27]. The FORCs, 
voltammetric currents, and FORC diagram are shown in Fig. [51 

Also indicated in Fig. [5|^a) are the FORC minima and the equilibrium 
isotherm, as calculated in independent equilibrium Monte Carlo simulations. 
Note that the FORC minima in Fig. [^(a) lie directly on the equilibrium 
isotherm. This is because such a system has one stable state for any given 
value of the potential, as defined by the continuous, single- valued equilibrium 
isotherm. The corresponding voltammetric currents are shown in Fig. [5](b). 
The uniformly positive value of the FORC diagram in Fig. M^c) reflects the 
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Figure 5: (Color online.) (a) FORCs for a continuous phase transition 
simulated at a slow scan rate, Q = 0.0003 meV/MCSS. The thin black curve 
near the middle of the major loop is the equilibrium isotherm. The inset is 
a magnification of the critical region. The minimum of each FORC is also 
shown (filled circles). The thick, black curve shows the first FORC which dips 
below the critical coverage, 9c ~ 0.36. (b) The corresponding voltammetric 
currents, calculated by numerical differentiation of the FORCs. See Eq. [T] 
(c) FORC diagram generated from the FORCs in (a). The positions of the 
FORC minima are shown as filled circles. The thick, straight line corresponds 
to the FORC marked as a thick curve in (a). 
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convergence of the family of FORCs with increasing /ij. This convergence 
results from relaxation toward the equilibrium isotherm, at a rate which in- 
creases with the distance from equilibrium. It is interesting to note that, 
while it is difficult to see at this slow scan rate, the rate of approach to 
equilibrium decreases greatly along the first FORC that dips below the crit- 
ical coverage 9^ ~ 0.36 (shown in bold in Fig. [S](a))- The FORCs that lie 
completely in the range 6 > 6c never enter into the disordered phase, and 
thus their approach to equilibrium is not hindered by jamming. This is a 
phenomenon that occurs when further adsorption in a disordered adlayer is 
hindered by the nearest-neighbor exclusion. As a result, extra diffusion steps 
are needed to make room for the new adsorbates, and the system follows 
different dynamics than a system with an ordered adlayer ^U\. The FORCs 
that dip below 6c enter into the disordered phase, and thus their approach to 
equilibrium is delayed by jamming. This is reflected in the FORC diagram 
by the Florida-shaped "peninsula" centered around this FORC in Fig. [5]|^c). 

The effect of jamming is more pronounced at higher scan rates, or with 
a higher diffusion barrier, where the rate of adsorption is much faster than 
the rate of diffusion. The family of FORCs and FORC diagram at a higher 
scan rate, Q = 0.01 meV/MCSS, are shown in Fig. [6l and with a larger 
diffusion barrier in Fig. [71 In Fig. [6l two distinct groups of FORCs undergoing 
jammed and unjammed dynamics can be clearly seen. This is reflected in the 
FORC diagram as a splitting of the "peninsula" into two "islands" of locally 
maximal values of p. A similar effect is seen in Fig. [71 since also there the 
rate of adsorption is much faster than the rate of diffusion (larger diffusion 
barrier). In addition. Fig. [Tl^a) shows a slight difference between the FORC 
minima and the equilibrium curve around the critical coverage. Notice also 
in Fig. M^a) that even at a much higher scan rate than in Fig. \5\ (nearly two 
orders of magnitude), the FORC minima still follow the equilibrium curve 
very accurately. Thus, the EC-FORC method should be useful to obtain the 
equilibrium adsorption isotherm quite accurately in experimental systems 
with slow equilibration rates. 

6 Comparison and conclusions 

Two observations can be made by comparing the FORCs, voltammetric cur- 
rents, and FORC diagrams for systems with discontinuous and continuous 
phase transitions. First, the FORC minima in systems with a continuous 
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Figure 6: (Color online.) (a) FORCs for a continuous phase transition 
simulated at a high scan rate, Q = 0.01 meV/MCSS, and the other model 
parameters as given in the text. The different curves and symbols have the 
same meanings as in Fig. [5](a). The inset showing the region of large fl and 
6 emphasizes the jamming behavior, (b) FORC diagram generated from the 
FORCs in (a). The different lines and symbols have the same meanings as 
in Fig. [5](b). 
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Figure 7: (Color online.) (a) FORCs for a continuous phase transition 
simulated with Q = 0.0003 meV/MCSS and a large diffusion barrier, Ann = 
300 meV. The other model parameters are as given in the text. The different 
curves and symbols have the same meanings as in Figs. [5](a) and[6](a). The 
inset showing the region of large /i and 6 emphasizes the jamming behavior, 
(b) FORC diagram generated from the FORCs shown in (a). The different 
lines and symbols have the same meanings as in Figs. EJ^b) and[6](b). 
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phase transition correspond closely to the equilibrium behavior, while they 
do not for systems with a discontinuous phase transition. Thus, FORCs 
can be used to recover the equilibrium behavior for systems with continu- 
ous phase transitions that need a long time to equilibrate. This should be 
useful in experiments. Second, due to the instability that exists in systems 
with a discontinuous phase transition, the minima of the family of FORCs in 
this case form a back-bending "van der Waals loop," and the corresponding 
FORC diagram contains negative regions that do not exist for systems with a 
continuous phase transition. Since experimental implementation of the EC- 
FORC method should only require simple reprogramming of a potentiostat 
designed to carry out a standard CV experiment, we believe the method can 
be of significant use in obtaining additional dynamic as well as equilibrium 
information from such experiments for systems that exhibit electrochemical 
adsorption with related phase transitions. 
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